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Projected-BCS wave functions have been proposed as the paradigm for 
the understanding of disordered spin states (spin liquids). Here we investi- 
gate the properties of these wave functions showing how Luttinger liquids, 
dimerized states, and gapped spin liquids may be described by the same 
class of wave functions, which, therefore, represent an extremely flexible 
variational tool. A close connection between spin liquids and "frozen" 
superconductors emerges from this investigation. 



Many years after the first proposal pj, the very existence of a spin- 
liquid ground state in two-dimensional (2D) spin-1/2 models is still a con- 
troversial issue. Short-range resonating-valence-bond (RVB) 2] phases 
with exponentially decaying spin-spin correlations and no broken lattice 
symmetry (i.e., with no valence-bond order), are conjectured to be stabi- 
lized by quantum fluctuations. However, while it is possible to show that 
spin-liquid ground states can be found in quantum dimer models Q, the 
numerical evidence in favor of disordered ground states in Heisenberg-like 
antiferromagnets is still preliminary 0, 0,1^, ■ Renewed interest on this 
topic is triggered by the possible realization of a frustrated Bose-Hubbard 
model with tunable interactions, by trapping cold atomic clouds in optical 
lattices 8] . A thorough analysis of the zero-temperature phase diagram of 
these models will probably require an outstanding numerical effort, how- 
ever, it is possible to tackle the spin-liquid problem by a different, less 
ambitious but remarkably informative, perspective. Here, we present a 
detailed analysis of the properties of a class of RVB wave functions, which 
have been shown to represent extremely good variational states for a wide 
class of microscopic Hamiltonians, ranging from effective low-energy mod- 
els for correlated electrons on the lattice to realistic models of atoms and 
molecules H El OH- In particular, we will show that both spin-liquid states 
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and valence-bond crystals may result from the "freezing" of Cooper pairs 
in the zero-doping limit of correlated electron models. 

Following Anderson's suggestion we define the class of projected- 
BCS (pBCS) wave functions in a TV-site spin lattice, starting from the 
ground state of a suitable BCS Hamiltonian: 

H(t, A) = ]T(i (J - fj,5 i:j ) c\ a Cj, a - [ A y c !,T c L + A u c ^ c *.T 

= Y^( 6k ~ & C L C M ~ [ A fc c l T c -fca + A fcC-fc,l C fe,T . (!) 

ka k 

where the bare electron band is real and both and are even 
functions of fc. A chemical potential fi is introduced in order to fix the 
number of electron equal to the number of sites. In order to obtain a class 
of non-magnetic, translationally invariant, singlet wave functions for spin- 
1/2 models, the ground state \BCS) of Hamiltonian is then restricted 
to the physical Hilbert space of singly-occupied sites by the Gutzwiller 
projector Pq. 

The first feature of such a wave function we want to discuss is the re- 
dundancy implied by the electronic representation of a spin state, i.e., the 
extra symmetries which appear when we write a spin state as Gutzwiller 
projection of a fermionic state. In turn, this property reflects in the pres- 
ence of a local, i.e., gauge, symmetry of the fermionic problem, as already 
pointed out several years ago j 1 1 1 ] . Let us consider a generic spin operator 
defined on a site in fermionic representation: X a p = c^cp. This operator 
acts in the Hilbert subspace of singly-occupied sites (which is left invari- 
ant under any spin Hamiltonian). It is easy to check that the three SU(2) 
generators N — (cjcf + cjcj, — l)/2, N + = cjct, A_ = (iV+)t commute, 
in the singly occupied site subspace, with X a p. This property reflects the 
invariance of the operator X a p under the usual U(l) gauge transformation 

94> '■ 

H ■ 4 - e»& (2) 
and also by the SU(2) rotation 

f cl\ ( cos 9 -isin6>\ f cl\ /„, 

These transformations are local, i.e., can be performed on each site in- 
dependently leaving every spin state invariant: they generate the SU(2) 
gauge symmetry group. 

Let us now consider the Gutzwiller projected BCS state in a Heisenberg- 
like model on a lattice with an even number N = In of sites. 

\ P BCS) = P G \BCS) = P G l[(u k +v k cl A cl kil )\0), (4) 

k 
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where the product is over all the N wave vectors in the Brillouin zone. 
The diagonalization of Hamiltonian Q gives explicitly 

I E k + e k A fc jEk — £fc / 2 i i a 12 

while the BCS pairing function is given by: 

, w fc A fc 

fk = — = — rrr- ( 5 ) 

Ufc £fc + iifc 

Clearly, the local gauge transformations previously defined change the BCS 
Hamiltonian, breaking in general the translation invariance. In the follow- 
ing, we will restrict to the class of transformations which preserve the 
translational symmetry of the lattice in the BCS Hamiltonian, i.e., the 
subgroup of global symmetries corresponding to site independent angles 
(0, 9). By applying the transformations J5J and iJS}, the BCS Hamiltonian 
keeps the same form with modified couplings: 

A^ -> A ije 2i * (6) 

for while the transformation Y,g gives: 

tij — > cos 20 — i sin cos 9 (Ay — A*, ) 

= cos 26> + sin 29 ImAy 

Ay -> (cos 2 6* Ay + sin 2 9 A*,- ) - i sin 29 t l} 

= RcA„ + i (cos 26* ImAy - sin 29 ) . (7) 

These relations are linear in tij and Ay and, therefore, equally hold for 
the Fourier components e k and Afc. Note that, being A r an even function, 
the real (imaginary) part of its Fourier transform Afc equals the Fourier 
transform of the real (imaginary) part of A r . It is easy to see that these two 
transformations generate the full 0(3) rotation group on the vector whose 
components are (efe, ReAfc, ImAfc). As a consequence, the length E k of this 
vector is conserved by the full group. In summary, this shows that there is 
an infinite number of different translationally invariant BCS Hamiltonians 
which, after projection, give the same spin state. Choosing a specific 
representation does not affect the physics of the state but changes the 
pairing function before projection. Within this class of states the only 
scalar under rotations, which can be given some physical meaning, is the 
BCS energy spectrum E k . Clearly, the projection operator will modify the 
excitation spectrum associated to the BCS wave function. Nevertheless, 
the invariance with respect to SU(2) transformations suggests that E k may 
reflect the nature of the physical excitation spectrum. 

Remarkably, in one dimension it is easy to prove that such a class of 
wave functions is able to faithfully represent both the physics of Luttinger 
liquids, appropriate for the nearest-neighbor Heisenberg model, and the 
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gapped spin-Peierls state, which is stabilized for sufficiently strong frus- 
tration. In fact, it is known [T^ that the simple choice of nearest-neighbor 
hopping tij (ek = —2t cosfc, fj, = ) and vanishing gap function A,j re- 
produces the exact solution of the Haldane-Shastry model, while choosing 
a next-nearest neighbor hopping (e/. = —2t cos2fc, /i = 0) and a sizable 
nearest-neighbor pairing (A = 4\/2f cos fc) we recover the Majumdar-Gosh 
state i.e., the exact ground state of the frustrated Heisenberg model 
when the next-nearest- neighbor exchange constant J2 is half of the nearest- 
neighbor one J\. Note that in this case, the BCS dispersion Ek is strictly 
positive, i.e., the BCS Hamiltonian is gapped. 

A quantitative analysis of this class of wave functions can be only carried 
out numerically. As an example, in Fig.l we show that the spin-spin cor- 
relations G(r) have a long-range tail, i.e., G(r) ~ 1/r, for both the known 
A = solution ^3] but also for our pBCS state with nearest-neighbor 
hopping and nearest- and third- neighbor A, whose energy is indeed re- 
markably accurate 0. Moreover, Fig. 2 shows that the same class of wave 
functions is able to display a clear dimer ordering as soon as the BCS dis- 
persion Ek shows a gap at the Fermi level. The agreement with the exact 
results (given by the Density-Matrix Renormalization Group method) for 
the frustrated Heisenberg model at J2/J1 = 0.4 is excellent also in this 
case. 
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FIG. 1: Spin-spin correlation function at the largest distance: G(r) = 
(— l) r (S*S*) in a chain with L=2r sites as a function of r. Circles: exact results 
for the Heisenberg model, Squares: pBCS wave function with nearest-neighbor 
hopping and third-neighbor A. Triangles: Gutzwiller projected Fermi sea. A 
1/r decay is clearly present in the variational wave functions. The exact result 
shows logarithmic corrections as predicted by conformal field theory. 

Now we specialize to the 2D square lattice and we investigate whether it 
is possible to further exploit the redundancy of the fermion representation 
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FIG. 2: Dimer-dimer correlation function (S^ +1 S^Sf So) as a function of a; cal- 
culated via the best variational pBCS function compared to the exact results ob- 
tained via Density-Matrix Renormalization Group (DMRG) in a 150-site chain. 

of a spin state in order to define a pairing function which, before projection, 
breaks the reflection symmetries of the lattice, while the projected state 
retains all the correct quantum numbers. We will show that, if suitable 
conditions are satisfied, a fully symmetric projected BCS state is obtained 
from a BCS Hamiltonian with fewer symmetries than the original spin 
problem. 

We first introduce a set of unitary operators. 

• Spatial symmetries: R x {x,y) — (x,—y); R xy (x,y) — (y,x). We 
define the transformation law of creation operators i?c] ff i? _1 = cjy 
and the action of the operator on the vacuum R |0) = |0). Note that 
these operators map each sublattice into itself. 

• Particle-hole: PhC^ a P h ^ 1 — i (— lpc^-o-, while the action of the Ph 
operator on the vacuum state is Ph\0) = Ylj ct ^.ct JO). 

• Gauge transformation: Get a G^ 1 = icj a and G |0) = |0). 

Clearly, R x and R xy are symmetries of the physical problem (e.g., the 
Heisenberg model) . G is a symmetry because the physical Hamiltonian has 
a definite number of electrons and Ph leaves invariant every configuration 
where each site is singly occupied if the total magnetization vanishes: N± = 
iVj = n. Therefore, the previously defined operators commute with the 
Heisenberg Hamiltonian and commute with each other because reflections 
do not interchange the two sublattices. The ground state of the Heisenberg 
model on a finite lattice, if it is unique, must be simultaneous eigenstate of 
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all the symmetry operators. We want to investigate sufficient conditions 
which guarantee that the projected BCS state is indeed eigenstate of all 
these symmetries. 

Let us consider a hopping term which just couples the two sublattices: 
e k+Q = —£k and a gap function A = A s +A a .2_ 2/ 2 +A xy with contributions 
from different symmetries. Moreover, we consider the case in which A s and 
A x 2_ y 2 will couple different sublattices while A xy is restricted to the same 
sublattice. If this is the case, the BCS Hamiltonian H(t, A s , A x z_ y 2, A xy ) 
transforms in the following way under the different unitary operators: 



R x H(t,A s ,A x 2_ y 2 7 A X y)R x 1 = H(t,A s ,A x 2_ y 2,—A xy ) 

R xy H(t,A s ,A x 2_ y 2,A xy )R X y — H(t 7 A s ,—A x 2_ y 2,A X y) 

P h H{t,A s ,A x 2_ y 2,A xy )P^ = H(t,A*,A x2 _ y2 ,-A* xy ) 

GH(t,A s ,A x 2^y2,A xy )G- 1 = H(t,-A s ,-A x 2_y2,-A xy ). 

Starting from these transformations, it is easy to define suitable compos- 
ite symmetry operators which indeed leave the BCS Hamiltonian invariant. 
For instance, let us consider the case in which A is real: in this case we 
select the R x Ph and R xy if A x 2_ y 2 = or R xy PhG if A s = 0. We cannot 
set simultaneously A x 2_ y 2 and A s different form zero. The eigenstates 
\BCS) of Eq. will be generally simultaneous eigenstates of these two 
composite symmetry operators with given quantum numbers, say a x and 
a xy . Let us consider the effect of projection over these states: 

a x P G \BCS) = P G R x P h \BCS)=R x P h P G \BCS) = 

= R X P G \BCS), (8) 

where we have used that both R x and Ph commute with the projector and 
that Ph acts as the identity on singly occupied states. Analogously, when 
a d x 2_ y 2 gap is present, 

a xy P G \BCS) = P G R xy P h G\BCS)=R xy P h GP G \BCS) = 

= {-l) n R xy P G \BCS). (9) 

These equations show that the projected BCS state with both xy and 
x 1 — y 2 contributions to the gap has definite symmetry under reflections, 
besides being translationally invariant. The corresponding eigenvalues, for 
n = N/2 even, coincide with the eigenvalues of the modified symmetry op- 
erators R x Ph and R xy PhG on the pure BCS state. Note that the quantum 
numbers we have defined, refer to the fermionic representation. It turns 
out that an extremely good variational wave function for the frustrated 
two dimensional Heisenberg antiferromagnet can be obtained by including 
gap functions of different symmetry 

An alternative, but equally interesting representation of the pBCS state 
can be given in terms of Slater determinants through the following argu- 
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ment. The pBCS wave function can be written in real space as: 



\ P BCS) = P G 



C R,a C X,- 



LR<X a 



|0) (10) 



where R and X run over the lattice sites and a is the spin index. The 
antisymmetric pairing function f(a,R;—a,X) = —f(—cr,X;a,R) is sim- 
ply related to the Fourier transform f(r) of the previously defined of 
Eq. f(^,R; l,X) = f(R— X). In order to avoid double counting, a 
given, arbitrary, ordering of the lattice sites has been assumed in Eq. (|10|> . 
By expanding Eq. i|l(J|) and defining f(a,R;a,X) = we get 

\ P BCS) = /{truRn^X!) 

Rl, (J\ ■ ■ ■ R n , On 

Xi,cf[ ■ ■ ■ X n , u n 

\oi,Ri;o[,Xi 

where the 2n labels (Ri ■ ■ ■ R n \ X\ ■ ■ ■ X n ) define a generic partition of the 
N = 2n sites of the lattice. We now rearrange the creation operators in 
the many body state of Eq. (| 1 1 1) 

\<j 1i R 1 \<t' 1i X 1 ■■■a n ,R n ;a' n ,X n ) = cj^^ct^, • ■ -ck^c^jO), 

according to a given site ordering in the lattice (irrespective of the spin) . 
This operation gives the (S z = 0) spin state \a\ ■ ■ ■ er/v) multiplied by a 
phase factor equal to the sign e of the associated permutation. However, 
many distinct terms correspond to the same spin state because they differ 
only by the way pairs are coupled: 



■ ■ ■ f(<j n ,R n ;a' n ,X n ) x 
• • • <7„, R„; a' , X n ) (11) 



\pBCS) = ^2^2e P f(a 1 ,R 1 ;a[,X 1 )---f(a n ,R n ;a' n ,X n ) | (Ti ■ • • a N ) 

Here P runs over all the possible (N — 1)!! partitions of the N = 2n sites of 
the lattices into pairs, for a given spin state |cri - - • <jn) uniquely identified 
by the variables {<7i}. The weight of the spin configuration is exactly the 
Pfaffian of the 2n x 2n antisymmetric matrix 14] : 



f(^R a -A,Rp) 



/(t, R a ] I, Xp) 
f(l,X a ;l,X p ) 



(12) 



where the matrix has been written in terms ofnxn blocks and R a are the 
positions of the up spins in the \<J\ - ■ ■ <Jn) state while X a are the positions 
of the down spins. The known relation 



[PfA] 



detA, 



(13) 
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together with the fact that we considered f(a,R;a,X) = imply that 

PfA = det f(R a — Xp) 0. Therefore, the weight of the spin state 

|<7i • • • a at) may be written as the determinant of a matrix which depends 
on the BCS pairing function. 

When we want to represent a spin state in a 2D lattice by the pBCS wave 
function, special care must be given to the definition of the phases of the 
pairing function f(r). As an example, let us consider a nearest-neighbor 
pairing function, i.e., we choose |/(f)| = 1 when r connects nearest neigh- 
bor sites in the lattice and zero otherwise. From Eq. Hl()[) it is clear that 
the pBCS state is written as the sum of all possible partitions of the N- 
site lattice into singlets (Ri,Xi) and the amplitude of a given partition is 
provided by the generic term of the Pfaffian of the matrix A: 

e P f(R 1 -X 1 )---f(R n -X n ). (14) 

Here we introduced the usual convention to orient a singlet from sublat- 
tice A to sublattice B and we fix Ri S A (Xi S B), while the permuta- 
tion P relates the chosen ordering of sites in the lattice and the sequence 
(Ri,Xi ■■ -R n ,X n ). Remarkably, as proved by Kasteleyn [14j. in planar 
lattices it is possible to choose the phase of the function f(r) so that the 
products of the form Eq. Ijl4|l have all the same sign. In such a case, 
the pBCS wave function exactly reproduces the short-range RVB state: 
the equal-amplitude superposition of all possible partitions of the lattice 
into nearest-neighbor singlets 0- In particular, on a rectangular lattice, 
the resulting pairing function has s + id symmetry fig, i.e., f(r) = 1 
on horizontal bonds and f(r) — i on vertical bonds. In turn, this may 
be obtained from the BCS Hamiltonian by considering third neighbor 
hopping ek = —2t (cos2k x + cos2fc y ) (with fj, = 0) and complex nearest- 
neighbor gap function = 8t (cos k x + i cos k y ). 

The properties of of the pBCS wave functions have been investigated by 
Lanczos technique and variational Monte Carlo method in the frustrated 
Heisenberg model on a square lattice, defined by the Hamiltonian: 

h = Ji J2 Si- Sj + J2Y. s > • s ;> ( 15 ) 

where the first sum runs over the nearest neighbors and the second on 
the next-nearest neighbors (i.e., along the diagonal). When both coupling 
constants are positive the model is frustrated. In Fig. 3 we compare the 
previously introduced dimer-dimer correlations (Si- Sf S^Sq) (where 
and (0, 1) are nearest neighbor sites) obtained via the optimized pBCS 
wave function and the exact Lanczos results for the Hamiltonian 1151) at 
J2/J1 = 0.55. The figure shows that the pBCS wave function is able 
to capture the correct behavior of correlations also in 2D. In Fig. 4 we 
report the size scaling of the squared dimer order parameter obtained via 
the optimal variational wave function in the pBCS class. The absence 
of dimer order is clearly suggested by the variational approach, pointing 
toward the existence of a spin-liquid in the 2D J1—J2 model. 
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FIG. 3: Dimer-dimer correlations in the 6x6 square lattice for the Hamilto- 
nian 11511 at J2/J1 = 0.55. Stars: exact Lanczos results. Triangles: optimized 
pBCS wave function. 
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FIG. 4: Size scaling of the squared dimer order parameter as predicted by the 
pBCS wave function. 



As a final remark, we like to comment on the fate of the insulating state 
described by a pBCS wave function upon doping. When a limited number 
of holes are injected into the lattice, it is likely that the basic structure 
of the pBCS state is not affected by the presence of mobile charges, be- 
ing largely determined by the super-exchange interaction among spins. As 
long as the BCS pairs present in the insulating state remain well defined 
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even at low doping, the system is generally expected to display super- 
conducting properties. This possibility has been explicitly verified in the 
so-called t—J model 9], where the natural generalization of a pBCS wave 
function has been shown to give rise to off-diagonal long-range order. Real 
materials, like high-temperature superconductors, display a considerably 
richer physics and other effects may inhibit the actual realization of this 
scenario: charge-density waves may occur at special values of the doping 
or the gain in hole kinetic energy may induce a global change in the spin 
state causing the breaking of electron pairs. Only a detailed study of the 
microscopic Hamiltonian will discriminate among these and other possibil- 
ities, and no general statement can be drawn on the basis of our analysis. 
However, our study suggests that a spin liquid (or even a valence bond 
crystal) may be thought of as an insulating state adiabatically connected 
to a superconducting phase which directly originates as soon as electrons 
are removed: this is a remarkable manifestation of the effects of electron 
correlations present in a gapped spin state which goes beyond standard 
band theory and the Fermi Liquid approach. 
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